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Abstract 

We  consider  the  use  of  preconditioning  methods  to  accelerate  the  convergence  to  a  steady 
state  for  both  the  incompressible  and  compressible  fluid  dynamic  equations,  VVV  also  consider 
the  relation  between  them  for  both  the  continuous  problem  and  the  finite  difference  approxi¬ 
mation.  The  analysis  relies  on  the  inviscid  equations.  The  preconditioning  consists  of  a  matrix 
multiplying  the  time  derivatives.  Hence,  the  steady  state  of  the  preconditioned  system  is  the 
same  as  the  steady  state  of  the  original  system.  For  finite  difference  methods  the  preconditioning 
can  change  and  improve  the  steady  state  solutions.  An  application  to  flow  around  an  airfoil  is 
presented. 


•T1  .4ii  research  was  partially  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA 
contract  Nos,  NASl-18605  and  NASl-19480  while  the  first  and  third  authors  were  in  residence  at  the  Institute  for 
Computer  Applications  in  Science  and  Engineering  (  ICA.SE),  NASA  Langley  Research  Center,  Hampton.  Va.  23681 


1  Introduction 


Seymour  Parter  has  considered  preconditioning  methods  for  numerical  approximations  to  elliptic 
partial  differential  equations.  As  an  extension  of  his  ideas  we  shall  consider  similar  techniques  for  the 
fluid  dynamic  equations.  Much  effort  has  been  expended  to  solve  the  compressible  steady  state  fluid 
equations  for  a  large  range  of  Mach  numbers.  A  standard  way  of  solving  the  steady  state  equations 
is  to  march  the  time  dependent  equations  until  a  steady  state  is  reached.  Since  the  transient  is 
not  of  any  interest  one  can  use  acceleration  techniques  which  might  destroy  the  time  accuracy  but 
enables  one  to  reach  the  steady  state  faster.  Such  methods  can  be  considered  as  preconditionings 
to  accelerate  the  convergence  to  a  steady  state.  For  the  incompresible  equations  the  continuity 
equation  does  not  contain  any  time  derivatives.  To  overcome  this  difficulty,  Chorin  [2]  added  an 
artificial  time  derivative  of  the  pressure  to  the  continuity  equation  together  with  a  multiplicative 
variable,  .  With  this  artificial  term  the  resultant  scheme  is  a  symmetric  hyperbolic  system  for  the 
inviscid  terms.  Thus,  the  system  is  well  posed  and  and  numerical  method  for  hyperbolic  systems 
can  be  used  to  advance  this  system  in  time..  The  free  parameter  (3  is  then  chosen  to  reach  the  steady 
state  quickly.  Later  Turkel  ([15],  [16])  extended  this  concept  by  adding  a  pressure  time  derivative 
also  to  the  momentum  equations.  The  resulting  system  after  preconditioning  is  no  longer  symmetric 
but  can  be  symmetrized  by  a  change  of  variables. 

It  is  well  known,  that  it  is  difficult  to  solve  the  compressible  equptions  for  low  Mach  numbers. 
For  an  explicit  scheme  this  is  easily  seen  by  inspecting  the  time  steps.  For  stability,  the  time  step 
must  be  chosen  inversely  proportional  to  the  largest  eigenvalue  of  the  system  which,  for  slow  flows, 
is  approximately  the  speed  of  sound,  c.  However,  other  waves  are  convected  at  the  fluid  speed,  u  , 
which  is  much  slower.  Hence,  these  waves  don’t  change  very  much  over  a  time  step.  Thus,  thousands 
of  time  steps  may  be  required  to  reach  a  steady  state.  Should  one  try  a  multigrid  acceleration  one 
finds  that  the  same  disparity  in  wave  speeds  s.uws  down  the  multigrid  acceleration.  With  an  implicit 
method  an  ADI  factorization  is  usually  used  so  that  one  can  easily  invert  the  implicit  factors.  The 
use  of  ADI  introduces  factorization  errors  which  again  slows  down  the  convergence  rate  when  there 
are  wave  speeds  of  very  different  magnitudes  [12]  . 

For  small  Mach  numbers  it  can  be  shown  ([5]>[7])[9])  that  the  incompressible  equations  ap¬ 
proximate  the  compressible  equations.  Hence,  one  needs  to  justify  the  computational  use  of  the 
compressible  equations  for  low  Mach  flows.  We  present  several  reasons  why  one  would  still  use  the 
compressible  equations  even  though  the  Mach  number  of  the  flow  is  small. 

•  There  are  many  highly  efficient  compressible  codes  available  that  could  be  used  for  such 
problems  especially  in  complicated  geometries. 

•  For  low  speed  aerodynamic  problems  at  high  angle  of  attack  most  of  the  of  the  flow  consists 
of  a  low  Mach  number  flow.  However,  there  are  localized  regions  containing  shocks. 

•  In  many  problems  thermal  effects  are  important  and  the  energy  equation  is  coupled  to  the 
other  equations.  Then,  the  compressible  equations  must  be  used  even  for  low  Mach  number 
flows. 

Therefore,  one  wants  to  change  the  transient  nature  of  the  system  to  remove  this  disparity  of  the 
wave  speeds.  Based  on  an  analogy  with  conjugate  gradient  methods  such  methods  were  called  [15] 
preconditioned  methods  since  the  object  is  to  reduce  the  condition  number  of  the  matrix.  Another 
approach,  in  one  dimension,  is  to  diagonalize  the  matrix  of  the  inviscid  term.  One  can  then  use  a 
different  time  step  for  each  equation,  or  wave.  Upon  returning  to  the  original  variables  one  finds 
that  this  is  equivalent  to  multiplying  the  time  derivatives  by  a  matrix.  Hence,  this  same  approach 
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is  named  characteristic  time  stepping  in  (17).  In  multidiinensions  one  can  no  longer  coinplelely 
decouple  the  waves  and  so  the  characteristic  time  stepping  is  only  an  approximation. 

Thus,  for  both  the  incompressible  and  compressible  equations  we  will  consider  systems  of  the 
form 

Wt  fx  +  9y  =  0. 

This  system  is  written  in  conservation  form  though  for  some  applications  this  is  not  necessary.  Our 
analysis  will  be  based  on  the  linearized  equations  so  that  the  conservation  form  Joes  not  appear 
in  the  analysis  though  it  does  appear  in  the  final  numerical  approximation.  This  system  is  now 
replaced  by 

P~'wt  +  /r  +  Py  =  0, 

or  in  linearized  form 

+  Au'x  4-  Bii-y  ^  C,  { I ) 

with  A  and  B  constant  matrices. 

In  order  for  this  system  to  be  equivalent  to  the  original  system,  in  the  steady  state,  we  demand 
that  P”^  have  an  inverse.  This  only  need  be  true  in  the  flow  regime  under  consideration.  We  shall 
see  later  that  frequently  P  is  singular  at  stagnation  points  and  also  along  sonic  lines.  Thus,  we 
will  temporarily  consider  strictly  subsonic  flow  without  a  stagnation  point.  For  transonic  flow  it  is 
necessary  to  smooth  out  the  singularity  in  a  neighborhood  of  the  sonic  Une.  We  also  assume  that 
the  Jacobian  matrices  A  =  and  B  =  ^  are  simultaneously  syminetrizable.  In  terms  of  the 
‘symmetrizing’  variables  we  also  demand  that  P  be  positive  definite.  We  shall  show  later,  in  detail, 
that  it  does  not  matter  which  set  of  dependent  variables  are  used  to  develop  the  preconditioner. 
One  can  transform  between  any  two  sets  of  variables.  Popular  choices  are  two  out  of  density, 
pressure,  enthalpy,  entropy  or  temperature  in  addition  to  the  velocity  components.  Thus,  when 
we  are  finished  we  will  analyze  a  system  which  is  similar  to  (1),  where  the  matrices  A  and  B  are 
symmetric  and  P  is  both  symmetric  and  positive  definite.  Such  systems  are  known  as  symmetric 
hyperbolic  systems.  One  can  then  multiply  this  system  oy  w  and  integrate  by  parts  to  get  estimates 
for  the  integral  of  wf,  i.e.  energy  estimates.  These  estimates  can  then  be  used  to  show  that  the 
system  is  well  posed  (see  e.g.  [5]).  We  stress  that  if  P  is  not  positive  then  we  may  change  the 
physics  of  the  problem.  For  example,  if  P  =  — /  then  we  have  reversed  the  time  direction  and  must 
therefore  change  all  the  boundary  conditions.  Hence,  to  be  sure  that  the  system  is  well  posed  with 
the  original  type  of  boundary  conditions  we  shall  only  consider  the  symmetric  hyperbolic  system. 

With  these  assumptions  the  steady  state  solutions  of  the  two  systems  are  the  same.  Assuming 
the  steady  state  has  a  unique  solution,  it  does  not  matter  which  system  we  march  to  a  steady  state. 
We  shall  later  see  that  for  the  finite  difference  approximations  the  steady  state  solutions  are  not 
necessarily  the  same  and  usually  the  preconditioned  system  leads  to  a  better  behaved  steady  state. 

2  Incompressible  equations 

We  first  consider  the  incompressible  inviscid  equations  in  primitive  variables. 

«r  +  Uy  =  0 

«i  +  UUx  +  VUy  +  Pr  =  0 

Vt  +  UVx  +  Vtty  +  Py  =  0 
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We  consider  generalizations  of  Chorin’s  pseudo-compressibility  method  [‘2],  Using  the  precondition¬ 
ing  suggested  in  [15]  (with  a  =  1)  we  have 


+  Ua:  -I-  t’y  =  0 

u 

^Pt  +  Ut+  UU^  -f-  VUy  4-  Pr  =  0 

^Pt  +  Vt  -f-  UUi  +  VVy  +  Py  =  0 


(2) 


or  in  conservation  form 


U 


-^Pt  -I-  tir  -h  Uy  -  0 

Pi  +  Ut  +  (tt*  +  p)j.  +  {Ut))y  =  0 

Pt  +  Vi  +  {uv)x  -f  -f  p)y  =  0 


We  can  also  write  (2)  in  matrix  form  using 


►-1 


/ 1//92 

0 

o\ 

( 

0 

0] 

1 

0  K 

p  = 

-u 

1 

0 

\ 

0 

1/ 

\ 

0 

1  / 

i.e. 


1//32  0  0 

u//?2  i  0 

W//32  0  1 


Multiplying  by  P  we  rewrite  this  as 


p  \  /  0  0  1 

u  -f-  I  0  V  0 
V  I  \  I  0  u 


=  0 


(3) 


wt  +  PAwx  +  PBwy  =  0. 


We  also  define 


D  —  A  + B  —  1  ^  u>i ,  u>2  ^  1 

where  wj ,  wj  are  the  Fourier  transform  variables  in  the  x  and  y  directions  respectively.  The  speeds 
of  the  waves  are  now  governed  by  the  roots  of  det{XI  -  P/ia>i  -  PBwa)  =  0  or  equivalently 
det(XP~^  -  A(4)i  -  Bu/3)  =  0.  Let 


q  =  uufi  +  vwj. 


Then  the  eigenvalues  of  PD  are 


<fo  =  9  (4) 

=s  ±/3 
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and  so  the  ‘acoustic’  speed  is  isotropic. 

The  spatial  derivatives  involve  symmetric  matrices,  i.e.  D  is  a  symmetric  matrix  but  P  is  not 
symmetric.  Thus,  while  the  original  system  was  symmetric  hyperbolic  the  preconditioned  system 
is  no  longer  symmetric.  In  [15]  it  is  shown  that  as  long  as 

0^  >  (u^  + 

then  the  equations  can  be  symmetrized.  On  the  other  hand  the  eigenvalues  are  most  equabzed  if 
4-  v^)  [15].  Hence,  we  wish  to  choose  j3^  ilightly  larger  than  +  v^.  However,  numerous 
calculations  verify,  that  in  general,  a  constant  /}  is  the  best  for  the  convergence  rate.  The  rea.sons 
for  this  are  not  clear. 

We  wish  to  stress  that  0  has  the  dimensions  of  a  spe«'d.  Therefore,  0  cannot  be  a  univer.'ial 
constant.  There  are  papers  that  claim  that  0  —  I  or  (i  =  2.5  are  optimal.  Such  claims  cannot  be 
true  in  general.  It  is  simple  to  see  that  if  one  nondimensionalizes  the  equation  then  0  gels  divided 
by  a  reference  velocity.  Hence,  the  optimal  ‘constant’  0  depends  on  the  dimensionalization  of  the 
problem  and  in  particular  depends  on  the  inflow  conditions.  In  many  calculations  the  inflow  mass 
flux  is  equal  to  1  or  else  p  +  +  u^)/2  =  1.  Such  conditions  will  give  an  optimal  0  close  to  one. 

However,  if  one  chooses  the  incoming  mass  flux  as  ten  then  the  optimal  0  will  be  larger. 

We  next  define  the  Bernoulli  function 

H  =  p  +  («=*  +  i;2)/2. 

Bernoulli’s  theorem  states  that  for  steady  inviscid  flow  H  is  constant  along  streamlines.  We  now 
multiply  the  second  equation  of  (2)  by  u  and  the  third  equation  of  (2)  by  v  and  add  these  two 
equations.  ]S  0^  =  ,  the  result  is 


Ht  +  uHx  +  vHy  =  0.  (5) 

Thus,  by  altering  the  time  dependence  of  the  equations  we  have  constructed  a  new  equation  in 
which  H  is  convected  along  streamlines.  Furthermore,  if  H  is  a  uniform  constant  both  initially  and 
at  inflow  then  H  will  remain  constant  for  all  time.  On  the  numerical  level  this  will  usually  not  be 
true  because  of  the  introduction  of  an  artificial  viscosity  or  because  of  upwinding.  For  viscous  flow, 
(5)  is  replaced  by 

Ht  +  uHx  +  vHy  —  -^{uAu  +  vAu) 

We  note  that  these  relationships  for  H  follow  from  the  momentum  equations  and  do  not  depend 
on  the  form  of  the  continuity  equation.  Hence,  we  consider  the  following  generalization  of  (2) 

-^pi  +  aHt  -i-Ux-^-Vy 

au 

~Pt  +  +  UUx  +  VUy  +  Px 

av 

-j^Pt  +  +  nVx  +  VVy  +  Py 

where,  a  is  a  free  parameter.  The  eigenvalues  of  PI)  are  independent  of  the  parameter  a  and  are 
given  by  (4).  For  a  =  0,0  =  1  we  recover  our  original  scheme.  For  a  =  -1  the  time  derivative  of 
the  pressure  no  longer  appears  in  the  continuity  equation.  For  general  0  we  have 


=  0 
=  0 
=  0 
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p-1 


(a  -i-  1)  au  av 
au  0'^  0 

av  0  0^ 


♦ 


f  f 


'■  =  5 


—au 


-au  l  +  a-ss^ 


-av 


aaui- 


-av 

gatiu 

l  +  a-^ 


! 


where  d  —  1  +  a  —  aa-^-  and  we  require  that  d  >  0.  If  0^  =  and  0  =  1  then 


/  +  „2 

:: 

In  [16]  an  analogy  to  the  symmetric  preconditioning  of  van  Leer,  Lee  and  Roe  was  constructed  for 
the  incompressible  equations.  If  we  choose  a  =  l,a  =  1  we  get  this  preconditioning  of  van  Leer 
et.al.,  i.e.  P  is  symmetric  . 

These  examples  show  that  the  preconditioning  is  not  unique.  If  fact,  since  the  determinant  of  the 
transpose  of  a  matrix  is  equal  to  the  determinant  of  the  original  matrix  it  follows  that  the  transpose 
of  P  is  also  a  preconditioner  with  the  same  eigenvalues  for  the  preconditioned  system.  In  general, 
these  various  systems  will  have  similar  eigenvalues  but  different  eigenvectors  for  the  preconditioned 
system.  Numerous  calculations  show  that  the  system  given  by  P  in  (2)  is  more  robust  and  converges 
faster  than  that  with  the  transpose  preconditioner.  This  shows  that  it  is  not  sufficient  to  consider 
just  the  eigenvalues  but  that  the  eigenvectors  are  also  of  importance.  However,  even  when  P  is 
symmetric  PD  is  not  symmetric  and  so  the  eigenvectors  of  the  preconditioned  system  do  not  form 
an  orthogonal  basis. 

We  next  examine  some  general  form  that  the  preconditioner  can  have.  For  this  analysis  it  is 
easier  to  use  strearawise  coordinates  as  suggested  in  [17]  and  so  u  =  0.  Let  w,  be  some  normalization 
of  the  velocity  components,  then 


/  0 

0  \ 

( 

0 

0 

u.  \ 

« 

II 

u 

0  , 

B  = 

0 

0 

0 

\  0 

0 

\ 

ti. 

0 

0  } 

Then  the  "convective”  eigenvector  for  the  non-preconditioned  system  is 


—au 


-av 

auv 


1  + 


The  "acoustic”  eigenvectors  are  given  by 


— uu(i  +v(tjwiF+^«2 
2 

U,Wi 

u,u;2 


—UU/I  -v/('iwi  )*+4u? 
2 

U,W2 
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We  now  consider  preconditioners  of  the  form 


a 

6 

0 

\ 

c 

d 

0 

0 

0 

1 

/ 

{(i) 


Let  D  =  ojjA  ,  u>f  +W2  —  1.  We  want  the  eigenvalues  of  PD  to  be  u^j  u.  ±u.  This  gives 

us  three  relations  for  the  four  unknowns: 


a  = 


(b  +  c)it,  +  dtt  =  0 
u^d-  bcvil  — 

The  values  suggested  in  [15]  are  6  =  0,c  =  =  1  while  the  values  suggested  in  [17]  are 

b  =  c  =  — =  2  .  We  next  present  the  eigenvectors  of  PD  in  terms  of  the  elements  of  P.  We 
exclude  the  case  uj2  —  0,u;j  =  ±1  as  in  this  case  PD  has  a  double  eigenvalue  u  and  the  eigensysteni 
completely  changes.  Then  the  "convective”  eigenvector  is 


The  "acoustic”  eigenvectors  are  given  by 


~  —  buu>?  —  +  buOJl 

u,(a  +  buf)  -  u,(6  +  c)a7i 
(Uthwi  +  «)u;2 


^  -  buuf  +  ^  -  btuji 
«,(a  +  buf)  +  u,(6  +  c)wi 
(utbufj  —  u)uj2 


We  note  that  the  convective  eigenvector  is  the  same  as  before  the  preconditioning  for  the  choice 
6  =  0.  The  two  acoustic  eigenvectors  are  orthogonal  to  each  other  if  we  choose  6  =  0  and  = 

—  .  This  is  similar,  but  not  identical,  to  the  choice  suggested  in  [1.5].  There  is  no  way  to  make 

the  convective  eigenvector  normal  to  both  acoustic  eigenvectors  for  preconditioners  of  the  form  (6). 


3  Compressible  equations 

The  time  dependent  Euler  equations  can  be  written  as 


Pt  +  UPx  +  VPy  +  pa^{Ul  +  Vy) 

Px 

Ui  +  UUx  +  VUy  +  — 

P 


Vt  + 


UVx  +  VVy  + 


Py 


St  +  uSx  +  vSy 


0 

0 

0 

0 


(7) 


where  a  is  the  speed  of  sound  given  by 

The  form  of  this  system  is  unchanged  if  we  nondimensionalize  the  equations.  From  now  on 
we  shall  assume  that  u,v,p,p  are  nondimensional  quantities  where  the  dimensional  variables  are 
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nondimensionalized  by  u.,p,,p,  with  p,  =  p.uj.  Following  [5]  we  define  ^  ^  •  If  the  fluid  is 

isentropic  then 

P=-4 


Hence,  as  e  goes  to  zero  the  speed  of  sound,  a,  goes  to  infinity  and  so  the  first  equation  in  (7) 
reduces  to  ttj,  +  Uy  =  0. 

It  was  pointed  out  in  ([15],  [16])  that  these  equations  can  be  symmetrized  by  using  ^  as  the 

independent  variable  rather  than  dp  .  Hence,  we  define  a  new  variable  <f)  by  d4>  =  ^.  For  isentropic 
flow  both  p  and  a  are  functions  only  of  the  density  and  so  using  (8,  9)  this  can  be  integrated 

iji 

explicitly.  This  gives  4>  =  .  As  the  Mach  number  goes  to  zero  4>  tends  to  infinity  and  therefore, 

Gustafsson  and  Stoor  [5]  subtract  a  constant  and  define 

x=i  , 

^  p  *  -  1 

2  * 

This  amounts  to  specifying  the  constant  in  the  integration  of  dip  from  dp.  They  then  prove,  using 
energy  methods,  that  for  the  linearized  equations 

,  dpincompretsikle 

- di - 

Since  p  -*  1  and  using  the  definition  of  d<p  this  is  equivalent  to 

d^Pcompreatible  dpincompreatikle  •  (19) 

We  consider  preconditionings  that  are  a  generalization  of  (3) 


fr  0  0  0 
^10  0 
^010 
0  0  0  1 


0  I  I  dt 

I  \dS 


/  «  «  0  0  \  /  g  \ 

I  a  u  0  0  du 

0  0  li  0  dv 

^  \0  0  0  u  I  \  dS  I  ^ 


^  u  0  o  0  ^ 

0  u  0  0 

a  0  V  0 

\  0  0  0  U  ; 


The  nonpreconditioned  case  corresponds  to  =  a^,a  =  0.  Let 


q  =  uui  +  vuf2 


then  the  eigenvalues  of  PD  are  given  by 


do  =  q  {double) 


d±=l/2  (1  —  a  +  la})q  ±  \J{{\  ~  a  ^  0'^ /a'^)‘^q'^  +  4(1  -  q'^/a'^)0'^ 

If  we  consider  the  special  case  a  =  1  +  /3^/a^  we  find  that  the  ‘acoustic’  eigenvalues  are  given  by 

d±  =  \/(l  - 


Hence,  these  eigenvalues  are  isotropic  in  the  limit  of  M  going  to  zero.  However,  this  eigenvalue 
vanishes  at  the  sonic  line  and  so  the  matrix  is  singular.  In  general  if  we  demand  that  the  acoustie 
eigenvalues  be  isotropic  then  we  have  a  singularity  at  the  sonic  line  where  the  eigenvalues  cannol 
be  isotropic.  The  two  ways  out  of  this  difficulty  are  either  to  smooth  the  formulas  near  the  singular 
line  or  else  to  give  up  on  isotropy.  This  difficulty  is  not  a  property  of  the  preconditioning  ju.st 
presented  but  applies  equally  to  all  preconditioners. 

We  now  consider  the  system  (7)  in  conservation  foru). 

pt  +  (p«)r  +  (/>w)v  =  0 

(pu)t  +  {pu^  +  p)i  +  {puv)y  -  0 

ipv)t  +  (put))r  +  {pv^  +  p)y  =  0 

Et  +  {pHu)j:  +  {pHv)y  =  0 

where 

_  p 

E  -  — —  + - 

7-1  2 

£  +  p 

”  "  —  =  + 

Note  that  the  Bernoulli  function  H  is  not  identical  with  H  for  the  incompressible  equations.  How 
ever,  we  again  have  that  for  steady  Inviscid  flow  H  is  constant  along  stream  lines.  We  now  pre¬ 
condition  the  density  and  the  energy  equations  in  the  following  consistent  manner.  Let  0  be  any 
variable  we  choose.  Then  we  consider 


V't  +  (p«)i:  +  {pw)y  =  0 
(t/j/f),  +  {pH  «)x  +  {pHv)y  =  0 

Manipulating  these  equations  gives 


Ht  +  uHx  +  vHy  =  0 


i.e.  the  total  enthalpy,  H,  is  simply  convected  in  time  along  streamlines  as  we  obtained  for  H  in  the 
incompressible  case.  It  is  interesting  to  observe  that  in  the  incompressible  case  we  achieved  this 
by  preconditioning  only  the  momentum  equations  while  for  the  compressible  flow  we  achieve  this 
by  preconditioning  the  continuity  and  energy  equation.  Of  course,  for  isentropic  flow  the  '^nergy 
equation  is  not  independent  of  the  other  equations  and  the  result  is  not  surprising. 

For  the  finite  difference  equation  this  implies  that  the  artificial  viscosity  for  the  continuity 
equation  should  be  based  on  ip  and  for  the  energy  equation  on  tj’H .  If  we  choose  V’  =  p,  i  e.  no 
preconditioning  for  the  continuity  equation  then  we  have  the  same  artificial  viscosity  as  suggested 
in  [6]  but  with  a  different  variable  being  advanced  in  uiae.  If  wo  choose  ij>  =  p  then  both  the 
continuity  and  energy  equations  are  preconditioned. 

We  next  present  the  van  Leer- Lee- Roe  preconditioning  for  general  non-aligned  flow  in  ( ^ ,  du,  di\dS} 
variables  [17]. 


/ 

-^u/a 
-■^vla 
\  0 


—  ^ti/ a 

0 


uv 


-jsv/a 

r  T  I  J  A 


(fTT  +  1) 


0 


+  T 


0 

0 

1  / 
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/  Vl  -  M2,  M  <  1, 

^  \  -  I,  M  >  1; 

/  >/l  -  M^,  M  <  1, 

^  j  \/l  -  A/'2,  M  >  1. 

At  the  sonic  line  (5  =  i'  ^nd  r  =  0  and  the  preconditioning  matrix  becomes  singular.  This 
preconditioning  is  not  unique  even  if  one  only  considers  symmetric  preconditioners.  In  both  these 
examples  the  preconditioner  was  constructed  based  on  using  (p,  u,  n,  .9)  a.s  the  dependent  variables. 
The  reason  for  this  choice  is  that  the  matrices  are  symmetric  which  this  choice.  However,  if 
another  '>>  ice  of  variables  is  more  appropriate  tht  introduces  no  difficulties.  Thus,  for  example 
[1]  recommends  the  use  of  (p,u,v,T)  variables  for  the  Navier-Stokes  equations.  Given  two  sets  of 
dependent  variables  w  and  W  let  be  the  Jacobian  matrix  ,  Then,  we  have  dW  =  W^,dv’.  So 
we  can  go  between  any  sets  of  primitive  variables  or  between  primitive  variables  and  conservation 
variables.  In  particular  since  the  equations  are  solved  in  conservation  variables  we  have  several 
ways  of  going  from  the  primitive  variable  preconditioner  to  a  conservation  variable  preconditioner 
Thus,  the  choice  of  variables  used  in  constructing  the  preconditioner  is  dictated  by  mathematical 
or  physical  reasoning  and  then  the  preconditioner  can  be  transformed  to  any  other  set  of  variables. 

•  Construct  the  preconditioner  matrix  for  the  conservation  variables.  If  W  are  the  conservative 
variables  and  w  the  primitive  variables  then  P conservative  =  {kI'u/)~'Pprim.«tvf(Ik’u,),  Details 
of  the  matrix  Jacobians  between  various  sets  of  variables  are  given  in  the  appendix. 

•  We  calculate  the  residual  dW  in  conservative  variables.  We  then  transform  dW  to  dxv  as 
before.  Next  we  multiply  by  P  and  finally  transform  back  to  conservative  variables  dW  and 
update  the  solution.  This  is  algebraically  equivalent  to  the  first  option  but  requires  three 
matrix  multiplies  instead  of  one.  However,  it  offers  more  flexibility. 

•  Similar  to  the  previous  suggestion  we  calculate  the  residual  dW  and  transform  to  conservative 
variables  dw  and  the  multiply  by  P.  At  this  stage  we  update  the  primitive  variables  w.  W'e 
then  use  the  nonlinear  relations  to  construct  W  from  w.  This  approach  has  advantages  if  the 
boundary  conditions  are  given  in  terms  of  the  primitive  variables  (p  or  T)  and  so  they  can 
be  specified  exactly  and  not  approximately. 

If  the  residual  dW  is  kept  from  the  conservation  form  but  the  time  derivative  Wt  is  replaced 
by  the  time  derivative  of  other  variables,  Wt  this  is  linearly  equivalent  to  preconditioning  by 
the  matrix  P~'  = 

oW 

These  methods  are  all  equivalent  for  linear  systems  and  the  difference  between  them  is  mainly 
one  of  convenience.  However,  we  shall  next  see  that  for  the  difference  approximation  these  ap¬ 
proaches  are  not  equivalent. 

4  Diflfereiice  Equations 

Until  now  the  entire  analysis  has  been  based  on  the  partial  differentia!  equation.  We  now  make 
some  remarks  on  important  points  for  any  numerical  approximation  of  thi,s  system. 

•  For  an  upwind  difference  scheme  based  on  a  Riemann  solver  this  Riemann  solver  .should  be 
for  the  preconditioned  system  and  not  the  original  scheme.  In  [3]  plots  are  .shown  to  illustrate 
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the  greatly  improved  accuracy  for  low  Mach  nmiiher  flows  when  the  Riematin  solver  is  based 
on  the  preconditioning. 

•  For  central  difference  schemes  there  is  a  need  to  add  an  artificial  vis-osity.  Act  iirai  y  is 
improved  for  low  Mach  number  flow’s  if  the  preconditioner  is  applied  (ndy  to  the  physical 
convective  and  viscous  terms  but  not  to  the  artificial  viscosity.  Volpe  [19]  siiows  that  the 
accuracy  of  the  original  system  deterioraU's  as  the  Mach  number  is  reduced,  i  he  use  of  a 
matrix  artificial  dissipation  ([14])  should  l)e  based  on  the  preconditioned  eipiations  as  in  tlse 
upwind  difference  scheme.  Upwind  schemes  without  preconditioning  teml  to  have  diffiruit ies 
with  accuracy  for  low  Mach  flows  [3], 

Hence,  both  for  upwind  and  central  difference  schemes  the  Riemaun  solver  or  artificial  visc<ts 
ity  should  be  based  on  P“^jP/lj  and  not  j/1j.  i.e.  in  one  dimen.sion  .solve  h’^+P =  (|P,4|iiv)j 
.  For  a  «calar  artificial  viscosity  |P-4|  is  replaced  by  the  spectral  radius  of  P.\  or  etpiivalentlv 
the  time  step  associated  with  the  preconditioned  matrix.  I  his  is  equivalent  to  not  multiplying 
the  artificial  viscosity  by  P. 

•  For  a  central  difference  scheme  with  a  scalar  artificial  viscosity  the  artificial  viscosity  is  of 
the  form  of  a  high  order  difference  of  the  same  quantity  as  is  advanced  in  time.  Thu.s,  the 
continuity  equation  is  solved  for  the  density  and  so  the  artificial  viscosity  i.s  a  difference  of  liie 
density.  Similarly,  for  the  momentum  equations.  For,  the  energy  etpiatioii  one  can  base  the 
artificial  viscosity  on  the  energy.  Alternatively  it  can  be  based  on  the  total  eiithalpy  which 
guarantee.s,  for  inviscid  flow,  that  the  total  enthalpy  is  constant  in  the  steady  state  [fij.  When 
preconditioning  the  system  one  of  the  alternatives  de.scribed  above  was  to  replace  tlie  time 
derivative  of  the  conservative  quantities  with  the  time  derivative  of  other  variables.  This, 
implies  that  the  artificial  viscosity  should  also  be  changed.  Thus,  if  the  continuity  equation 
is  updated  for  the  pressure  rather  than  the  density,  then  tiie  artificial  viscosity  should  he 
based  on  the  pressure.  This  i.s  physically  more  reasonable  for  low  speed  flow  since  tin*  density 
is  almost  constant  and  so  will  not  contribute  any  reasonable  viscosity.  Furthermore,  using 
a  vi.scosity  in  the  continuity  equation  based  on  the  pressure  mimics  what  was  done  for  the 
incompressible  equations.  This  allows  the  low  speed  compressible  equations  to  replicate  the 
results  of  the  incompressible  equations  on  the  finite  difference  level.  This  will  be  (.isrussed  in 
more  detail  in  th«  following  .sections. 

•  When  using  characteristics  in  the  boundary  conditions  these  should  bp  based  on  the  ciiarac- 
teristics  of  the  modified  system  and  not  the  physical  system. 

•  When  using  multigrid  it  is  better  to  transfer  the  residuals  based  on  the  preconditioned  system 
to  the  next  grid  since  these  residuals  are  more  balanced  than  the  physical  residuals. 

Preconditioning  is  even  more  important  when  using  multigrid  than  with  an  explicit  .scheme. 
W^ith  the  original  system  the  disparity  of  the  eigenvalues  greatly  affects  the  smoothing  r;ites 
of  the  slow  components  and  .so  slows  down  the  multigrid  method,  [10]. 

•  As  indicated  above  there  are  accuracy  difficulies  at  low  Mach  numlaT.s  [19].  .Some  of  these 

can  be  alleviated  by  preconditioning  the  dissipation  terms.  For  very  small  Mach  numbers 
there  is  also  a  difficulty  with  roundoff  errors  as  —  oo.  Several  [leople  have  stigge.sied 

subtracting  out  a  constant  pressure  from  the  dynamic  pressure.  A  more  detailed  analysis 
[4]  sugge.sts  replacing  the  pressure  p  by  p  where  p  ~  and  c  is  a  representative  .Vlach 

number. 
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We  conclude  from  the  above  remarks  that  the  steady  state  solution  of  the  preconditiom-d  ^yh 
teni  may  be  different  from  that  of  the  physical  system.  Thus,  on  the  finite  difference  level  itie 
preconditioning  can  improve  the  accuracy  as  well  as  the  convergence  rate. 

In  the  previous  section  we  stated  that  it  is  not  important  if  one  updates  a  different  set  of  vari 
ables  or  else  uses  the  conservation  variables  and  compensates  with  preconditioning  by  a  matrix 
multiplication.  However,  numerically  for  very  small  Mach  numbers  the  eiitrir.s  in  the  precondiiiun 
ing  matrix  can  become  very  large  or  small.  Hence,  it  can  be  advantageous  to  update  the  pressure 
or  temperature  directly  rather  than  using  a  matrix  multiply  for  preconditioning. 

5  Convergence 

We  have  previously  quoted  several  papers  ((•'j],(7j.  [9j),  that  prove  the  convergence  of  the  compress 
ibie  equations  to  the  incompressible  equations,  for  i.sentropic  flow,  as  the  Mach  miiiiber  goes  to 
zero.  For  nonisentropic  flow  there  are  no  formal  proofs.  However,  it  is  clear  that  for  viscous  flows 
that  the  boundary  condition  on  the  temperature,  adiabatic  or  isothermal  is  very  important,  see 

[Ill- 

All  these  results  refer  to  the  time  dependent  physical  equations.  Once  preconditioning  is  intro¬ 
duced  time  accuracy  is  lost  and  one  can  only  discuss  convergence  of  the  steady  state  solutions.  In 
this  case  one  would  hope  that  the  time  dependent  preconditioned  compressible  equations  converge 
to  some  time  dependent  preconditioning  of  the  incompres.sible  equations.  In  addition,  one  would 
also  want  this  to  be  true  on  the  numerical  level.  Thus,  one  would  want  to  solve  the  preconditioned 
compressible  pquation.s  by  some  numerical  technique,  on  a  fixed  mesh  and  compare  that  with  the 
solution  of  the  incompressible  equations  on  the  same  mesh.  .Mathematically,  we  have  two  limit  pro¬ 
cesses  occuring:  the  Mach  number  going  to  zero  and  the  mesh  size  going  to  zero.  The-se  two  limits 
need  not  commute.  If  one  first  converges  the  mesh  size  and  then  the  Mach  number  it  is  equivalent 
to  the  convergence  proofs  for  the  analytic  case.  The  more  interesting  problem  is  to  converge  the 
Mach  number  and  then  converge  the  mesh,  i.e.  we  u.se  a  fix  me,sh  as  the  Mach  number  is  reduced. 

In  particular  this  requires  a  careful  study  of  the  viscosities  introduced  by  the  scheme.  We  first 
consider  an  upwinding  scheme.  For  the  compressible  case  we  have  already  noted  that  the  Rieniann 
solver  should  depend  on  the  preconditioned  problem.  Ojie  would  then  need  to  .show  that  this  Rie- 
mann  problem  converges  to  a  Riemann  problem  for  some  preconditioning  of  the  incompressible 
equations.  We  next  consider  a  central  difference  .scheme  with  a  .scalar  viscosity.  In  this  case  a  high 
order  even  difference  of  some  quantity  is  added  separately  to  eadi  equation,  e.g.  for  the  incoiti- 
pre.ssible  equations:  pressure  for  the  continuity  equation,  u  and  v  for  the  momentum  equations.  For 
the  compressible  equations  one  normally  adds  a  density  difference  to  the  continuity  equation.  In 
such  a  case  it  is  obvious  that  the  numericai  scheme  for  the  compressible  etpiatiniis  cannot  converge 
to  the  numerical  scheme  for  the  incompre.s.ssible  equations.  Furthermore,  for  low  .Mach  number 
flows  the  density  is  almost  constant  and  so  the  higher  order  difference  of  the  density  does  not  add 
much  of  a  viscosity  to  the  continuity  equation.  .As  such,  we  conclude  that  the  artificial  viscosity 
for  the  compressible  continuity  eq\iation  should  be  based  oti  pre.ssure  and  not  density  (at  least  for 
low  Mach  numbers). 

We  shall  examine  the  convergence  a  little  more  closely.  Hy  convergence  (>f  the  compressible 
equations  to  the  incompressible  eqtiations  we  are  merely  verifying  what  happens  to  the  difference 
equations  a.s  the  .Mach  number  goe.s  to  zero.  The  convergence  of  the  solution  of  the  numerical 
approximation  to  the  preconditionined  compressible  equations  to  the  nnmerical  solution  of  the 
incompressible  eqtiation  is  more  difficult.  However,  we  shall  see  that  for  the  numerical  solution  the 
convergence  of  the  difference  equation  i.s  nontrivial  and  depends  on  the  preconditioning.  For  this 
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purpose  we  shall  only  consider  a  central  difference  approximation  together  with  a  scalar  artificial 
viscosity  for  the  nondimensionalized  preconditioned  inviscid  equations. 

For  the  incompressible  equations  in  nonconservative  form  we  consider  the  preconditioned  system 


Pt  +  4-  Vy) 

u 

■^Pt  +  Uj  +  UUx  +  VUy  +  Px 
■^pt  +  Ut  +  UVx  VVy  +  Py 


^  I(  f^”lPrrr  )x  4"  (  A  2Pj(yy  )  j) 
Wx3;x)r  4-  (  A  itiyyj,  )j,] 

A  [(  A  1  nxrr)x  4"  {  A  2t'yyy  )y] 


(H) 


where  each  space  derivative  is  approximated  by  a  central  difference  with  spacing  h  in  each  direction. 
The  time  derivatives  are  replaced  by  a  multi-stage  scheme.  A'l,  K-z  are  the  largest  eigenvalues  of 
the  coefficient  matrix  in  the  respective  direction.  Since  we  do  not  expect  shocks  we  only  consider 
a  linear  fourth  difference  artificial  viscosity  and  not  a  nonlinear  second  difference  [6],  see  the  resuit 
section  for  more  details. 

We  next  consider  the  same  scheme  foi  tlie  preconauioned  compressible  inviscid  equations,  under 
the  assumption  that  the  entropy,  S,  is  constant  so  that  p  =  p(p).  It  easier  to  analyze  the  convergence 
for  the  nonsymmetric  form  since  the  pressure,  p,  convergences  and  not  see  (10).  For  the 
preconditioned  continuity  equation  we  have 

Pi  4-  ^  \upx  +  VPy  -h  pa}{Ux  4-  Vy)j  =  [(  A']PxXx)x  4-  (  A'2Pyyy)y] 

Since  p*  =  a^p^iPy  =  “^Pv  can  rewrite  the  system  as 

Pi  +  [{pu)x  4-  (pw)y))  =  ((A'lPrrx)x  4-  iK-2Pyyy)y] 

li  "px 

~Pl  4-  U(  4-  UUx  +  VUy  +  ^  =  [iKiUx=x)r  +  ( A'2Myyv)y]  ( 

■J^Pt  +  Vt  UVx  +  VVy  +  ^  =  A^[(A'lV-Xl)x  4-  (A'2Uyyy)y] 

Comparing  (11)  with  (12)  it  is  obvious  that  if  p  — ►  1  as  M  — ►  0  then  (11)  converges  to  (12). 
It  is  crucial  for  both  the  time  derivative  and  the  artificial  viscosity  in  the  compressible  continuity 
equation  to  be  pressure  based  rather  than  density  based.  The  preconditioning  of  the  momentum 
equations  is  not  important  for  this  convergence. 

For  the  incompressible  equations  in  conservative  form  we  multiply  the  first  equation  in  (11)  by 
u  and  add  it  to  the  second  and  third  equations.  However,  we  do  not  change  the  artificial  viscosity. 
Then 


Pi  4-  /'  ’lx  -f-  Vy) 
~Pl  +  Ut  +  (u^  +  p)x  4  («1.’)y 
2  V 

IjyF*  4-  fi  4-  (uv)x  +  (v^  4-  p)y 


h'  [{  A  tPixx)x  4"  (  A  2Pyvi  )yj 

A  [(  A  ]  +  (  A  2Uyyy  )y] 

h'  [(  A  ]  I'rxx  )x  4"  (  A  2  Vyyy  )y  j 


(13) 


For  the  compressible  equations  in  conservative  form  we  have  two  choices.  One  choice  is  to 
multiply  the  first  equation  in  (11)  by  u  and  the  second  by  p  and  add  the  two.  The  spatial  derivatives 
are  then  in  conservation  form.  However,  the  time  derivative  is  of  the  form  p«(  rather  than  (pu); 
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and  the  artificiaj  viscosity  terms  are  not  in  conservation  form.  Hence,  we  instead  choose  to  apply 
the  preconditioning  directly  to  the  conservative  form.  The  resultant  preconditioned  compressible 
equations  in  conservative  form  is 


pt  +  0^  [(pw)*  +  (pn)s,)j 
2u 

—Pt  +  {fm)t  +  (pti^  +  p0  +  ipuv)y 
~Pt  +  {pv)t  +  (pnn)x  +  (pv^  +  p)y 


h  iPTxi:)z-  if^2Pyyy)y] 

h^Kiipu) 

rxx)r  "b  {  A  i{p'^)yyy)y\ 
h^[{K,{pv)  XXX  )x  +  (A'2(pn)yj,y)y] 


(14) 


Note  that  (14)  is  not  equivalent  to  (12). 

In  this  case  we  again  see  that  (14)  converges  formally  to  (13)  as  A/  -+  0  and  p  — »  1.  This  is 
because  the  pressure  is  used  for  the  time  derivative  and  the  artificial  viscosity  in  the  continuity 
equation. 

This  aU  applies  to  the  isentropic  equations.  The  compressible  equations  for  nonisentropic  flow 
is  more  complicated  and  in  fact  there  does  not  exist  any  proof  of  the  convergence  of  the  solution 
of  the  compressible  equations  to  the  solution  of  the  incompressible  equations  for  this  case. 


6  Computational  Results 

We  now  present  a  calculation  for  two  dimensional  flow  around  an  airfoil  to  demonstrate  the  previous 
theory.  As  described  above  the  discretization  is  based  on  the  multistage  time  method  coupled  with 
a  central  difference  approximation  as  described  in  [6]. 

We  solve  the  equation  in  conservation  form  based  on  a  hybrid  set  of  variables  of  those  previously 
considered. 


Wt  +  P(Fx  +  Gy)  =AD  =  (A',(?„x)r  +  {K2Qyyy)y 

W=\  I" 

pv 


<  pu  ^ 

''  pu  ^ 

(  P'  \ 

F  = 

pu^  +  0 

.  G  = 

puv 

,  Q  = 

pu 

puv 

pv^  +  p' 

pv 

\  pB'u  j 

1 

\  pB'v  ) 

[  H'  ) 

where  p'  =  p  —  p^,  E'  =  Cpp{T  -  T^o)  —  {p  -  Poo)  +  Ip’  )  and  pH'  =  E'  +  p'  .  We  subtract  the 
constants  to  keep  the  quantities  in  scale,  see  (10). 
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We  choose 


B2 

B3 


Ba 


Biii  + 


Biv  + 

B^ff  + 


au 

w 

av 


1''  '  02 

a(u^  4-  v^) 


/3*  =  max(u^  +  v®,0.9(«^  +  v^)),  a=l 

These  equations  are  given  for  for  the  nondimensionalized  variables.  The  nondiinensionalization 
affects  the  convergence,  In  some  codes,  p  and  p  are  fixed  in  the  far  field.  This  implies  that  the 
speed  of  sound,  a,  is  also  bounded.  As  the  Mach  number  goes  to  zero  the  pressure  remains  of  order 
1  while  the  velocities  go  to  zero.  Alternatively,  one  can  nondimensionalize  so  that  the  velocities 
are  of  order  1  in  the  far  field  and  then  the  pressure  and  speed  of  sound  go  to  infinity,  unless  one 
subtracts  an  appropriate  constant, 

A  typical  step  of  a  Runge-Kutta  approximation  is 

=  1^(0)  _  -  AD\  , 

where  Dx  and  Dy  are  spatial  differencing  operators,  and  AD  represents  the  artificial  dissipation 
terms.  The  dissipation  terms  are  a  blending  of  second  and  fourth  differences.  That  is, 

AD={Dl-^Dl-Dt-D^^Q,  (15) 


where 


DiQ = V.  ^-]  «'■'' 

DiQ  =  V. 


and  Aj;,Vi  are  the  standard  forward  and  backward  difference  operators  respectively  associated 
with  the  X  direction.  The  variable  scaling  factor  A  is  chosen  as 


-  2 


where  A*  and  Xy  are  proportional  to  the  largest  eigenvalues  of  the  matrices  A  and  B.  For  genei  alized 
coordinates  x  and  y  are  replaced  by  ^,77  respectively.  This  spectral  radius  is  now  a  function  of  the 
preconditioning.  Hence, 

A,.  =  p(PA)  A„  =  p(PB) 

The  coefllcients  and  are  adapted  to  the  flow  and  are  defined  as  follows: 


=  max(t/,_i  j,  I'ij,  t't+2,i), 


= 


Pt+l,j'  ~  ^Pi,j  ~h  Pi-1,; 
Pi+l,i  +  2pi,j  +  Pi-1, j 


max 
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where  p  is  the  pressure,  and  the  quantities  and  are  constants  to  be  specified.  The  operators 
in  (15)  for  the  y  direction  are  defined  in  a  similar  manner. 

The  second-difference  dissipation  term  is  nonlinear.  Its  purpose  is  to  introduce  an  entropy-like 
condition  and  to  suppress  oscillations  in  the  neighborhood  of  shocks.  This  term  is  small  in  the 
smooth  portion  of  the  flow  field.  The  fourth -difference  dissipation  term  is  basically  linear  and  is 
included  to  damp  high-frequency  modes  and  aUow  the  scheme  to  approach  a  steady  state.  Ordy  this 
term  affects  the  linear  stability  of  the  scheme.  Near  shocks  it  is  reduced  to  zero.  For  incompres.siljh* 
flow  shocks  can  only  appear  in  the,  nonphysical,  transient  and  so  the  second-difference  dissipation 
is  not  important.  To  reemphasize,  the  preconditioning  matrix  multiplies  the  the  flux  terms  but  not 
the  artificial  viscosity  terms.  The  scaling  in  the  artificial  viscosity  depends  on  the  spectral  radius  of 
the  preconditioned  matrices.  If  one  were  to  use  a  matrix  valued  viscosity,  [14]  ,it  would  be  related 
to  the  absolute  value  of  the  preconditioned  Jacobian  matrices. 

The  boundary  conditions  at  the  far  field  boundary,  for  subsonic  flow,  are  based  on  the  one 
dimensional  theory  of  characteristics  in  the  direction  normal  to  the  boundary.  The  preconditioning 
changes  the  form  of  these  characteristic  variables.  They  are  now  given  by 


where  u  is  the  component  of  the  velocity  normal  to  the  boundary.  This  formulas  simplify  slightly 
if  a  =  1  and  more  if  a  =  1  -I-  |r-  If  we  consider  low  Mach  numbers  then  we  can  approximate  these 

by 


=  U - a’ 


R2 


U  -f 


p(i 


which  is  the  same  as  for  the  incompressible  case.  At  solid  boundaries  the  normal  momentum 
equation  is  used  which  is  not  affected  by  the  preconditioning. 

The  solution  is  advanced  by  the  explicit  Runge-Kutta  method  described  above  and  without 
any  residual  smoothing  or  multigrid.  We  present  two  calculations  for  inviscid  flow  about  a  NACA 
0012.  The  first  calculation  is  for  inflow  conditions  M  =  0.03,  a  =  4°  .  In  this  case  we  see 
that  the  residual  asymptotes  without  the  use  of  preconditioning  and  that  the  preconditioning 
dramatically  increases  the  rate  of  convergence.  The  use  of  the  preconditioning  adds  only  a  few 
percent  to  the  total  computational  time.  For  viscous  flows  the  computational  time  required  for 
the  preconditioning  is  negligible.  In  the  second  case  we  consider  the  same  geometry  but  with  an 
inflow  of  M  =  0.8,  a  =  1.25®  .  In  this  case  we  also  see  a  increased  rate  of  convergence  for  the 
preconditioned  case  but  not  as  dramatic  as  before. 

In  all  cases  we  could  not  allow  /?  to  become  too  small.  In  fact  the  cutoff  is  sufficiently  large 
so  that  ()  is  close  to  a  constant.  This  has  been  observed  in  many  central  difference  Runge-Kutta 
codes  but  has  not  been  observed  in  the  upwind  code  coupled  with  an  ADI  solver  [3]  . 


7  Conclusion 

We  have  considered  a  family  of  matrix  preconditionings  for  both  the  incompressible  and  compress¬ 
ible  fluid  dynamic  equations  that  generalize  previous  results.  In  both  cases  the  wave  speeds  are 
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more  equalized  than  for  the  original  set  of  equations  and  so  the  condition  number  of  the  system  is 
reduced.  For  the  compressible  equations  the  condition  is  equal  to  1  at  a  Mach  number  of  zero  and 
increases  as  the  Mach  number  increases.  At  M  =  1  the  condition  number  is  infinite  but  it  increases 
at  a  slower  rate  than  for  the  physical  system. 

In  addition  to  the  question  of  the  convergence  rate  to  a  steady  state  we  have  considered  the 
question  of  the  accuracy  of  the  numerical  scheme  for  low  Mach  numbers.  One  can  prove  that  for  the 
partial  differential  equation  that  the  compressible  equations  approach  the  incompressible  equations 
as  the  Mach  number  goes  to  zero.  For  the  numerical  scheme  this  is  no  longer  generally  true  and  so 
the  accuracy  of  the  numerical  scheme  to  the  compressible  equations  decreases  as  the  Mach  number 
goes  to  zero.  One  way  to  improve  the  situation  is  to  include  the  preconditioning  in  the  Riemann 
solver,  or  equivalently,  to  account  for  the  preconditioning  in  the  artificial  viscosity.  For  example, 
for  low  Mach  numbers  the  scalar  artificial  viscosity  for  the  continuity  equation  should  be  based  on 
the  pressure  rather  than  the  density. 
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A  Appendix 


Let  W  denote  the  conservative  variables  (p,  m,  n,  £)‘,  with  m  =  pn,n  =  pv  ,  let  w  denote  the 
primitive  variables  (p,  w,v,5)‘  and  let  w  denote  (p,  u,n,T)L  Then 


dw 

dw 


( 


\ 


j 

ti 

?■ 

M 

?■ 

1  ,  M* 

-r-l  2 


0  0 
9  0 
0  p 

m  n 


Cp 

_  S. 

_ ^ 

Cp  I 

2cp  / 


dw 

w 


2 

_  3i 

P 

_ v 

<=p(('r-i)V-i] 
\  - - 


-(7  -  1)“ 

i 

o 

0 


7-1  '' 

0 
0 


/ 


1 

e. 

p 

0 

0 

dw  _ 

m 

p 

9 

0 

m 

~T 

dw 

n 

p 

0 

9 

n 

~T 

\ 

s. 

p 

m 

n 

p(u^  -f-v*  ) 

2T 

\ 


/ 


2 

-(7-l)u  -(7-l)v 

7-l\ 

dw 

—  a 

p 

0 

dw  ~ 

—  Ji 

p 

0  i 

p 

0 

r  , 

\  p  ^  2p 

(-11-1)11  (7-l)v 

pR  pR 

/ 

References 

[1]  Choi,  Y.-H.  and  Merkle,  C.L.,  The  Application  of  Preconditioning  to  Yiscous  Flows,  to  appear 
JournaJ  of  Computational  Physics. 

[2]  Chorin,  A.J.,  v4  Numerical  Method  for  Solving  Incompressible  Viscous  Flow  Problems,  Journal 
of  Computational  Physics,  2  (1967),  pp.  12-26. 

[3]  Godfrey,  A.G.,  Walters.  R.W.  and  Van  Leer,  B.,  Preconditioning  for  the  Navier-Stokcs  Equa¬ 
tions  with  Finite-Rate  Chemistry,  AIAA  paper  93-0.53.'),  1993. 

[4]  Gustafsson,  B.,  Unsymmetric  Hyperbolic  Systems  and  the  Euler  equations  at  Low  Mach  Num¬ 
bers  ,  J.  Scientific  Computing,  2(1987),  pp.  123-1.36. 

[.5]  Gustafsson,  B.,  and  Stoor,  H.,  Navier-Stokcs  Equations  for  Almost  Incompressible  Flow,  SIAM 
J.  Numer.  Anal.,  28(1991),  pp.  1.523-1547. 

[6]  Janie.son,  A.,  Schmidt,  W.  and  Turkel,  E.,  Numerical  Solutions  of  the  Euler  Equations  by 
a  Finite  Volume  Method  using  Runge-Kutta  Time-Stepping  Schemes  ,  ALA  A  paper  .^1-1259, 
1981. 


17 


[7]  Klainerman,  S.  and  Majda,  A.,  Compressible  and  Incompressible  Fluids  ,  Comm.  Pure  Appl. 
Math.,  35(1982),  pp.  629-651. 

[8]  Kreiss,  H.O.  and  Lorenz,  J.,  Initial-Boundary  Value  Problems  and  the  Navier-Stokes 
Equations,  Academic  Press,  Boston  ,  1989. 

[9]  Kreiss,  H.O.,  Lorenz, J.  and  Naughton,  M.,  Convergence  of  the  Solutions  of  the  Compressible 
to  the  Solutions  of  the  Incompressible  Navier-Stokes  Equations  ,  Advances  in  Appl.  Math., 
12(1991),  pp.  187-214. 

[10]  Len,  D.  and  Van  Leer,  B.,  Progress  in  Local  Preconditioning  of  the  Euler  and  Navier-Stokes 
Equations,  11th  AIAA  CFD  Conference,  AIAA  paper  93-3328,  pp.  338-3481993. 

[11]  Panton,  R.L.,  Incompressible  Flow,  pp.  237-260  ,  Wiley-Interscience,  John  Wiley  and  Sons, 
Neiv  York,  1984. 

[12]  Steger,  J.L.  and  Kutler,  P.,  Implicit  Finite- Difference  Procedures  for  the  Computation  of  vortex 
wakes  ,  AIAA  J.,  15(1977),  pp.  581-590. 

[13]  Swanson,  R.C.,  Turkel,  E.,  Pseudo-Time  Algorithms  for  Navier-Stokes  Equations  ,  Applied 
Numerical  Math.,  2(1986),  pp.  321-334. 

[14]  Swanson,  R.C.,  Turkel,  E.,  On  Central  Difference  and  Upwind  Schemes  ,  Journal  of  Compu¬ 
tational  Physics,  101(1992),  pp.  292-306. 

[15]  Turkel,  E.,  Preconditioned  Methods  for  Solving  the  Incompressible  and  Low  Speed  Compressible 
Equations  ,  Journal  of  Computational  Physics,  72(1987),  277-298. 

[16]  Turkel,  E.,  A  Review  of  Preconditioning  Methods  for  Fluid  Dynamics,  to  appear  Applied 
Numerical  Mathematics. 

[17]  Van  Leer,  B.,  Lee,  W.T.,  Roe,  P.L.  ,  Characteristic  Time-Stepping  or  Local  Preconditioning  of 
the  Euler  Equations  ,  AIAA  Paper  91-1552, 1991. 

[18]  Van  Leer,  B.,  Lee,  W.T.,  Roe,  P.L.,  Powell,  K.G.,  Tai,C.-H.,  Design  of  Optimally  Smoothing 
Multistage  Schemes  for  the  Euler  Equations ,  Comm.  Applied  Numer.  Meth.,  8(761-769),  1992. 

[19]  Volpe,  G.,  Performance  of  Compressible  Flow  Codes  at  Low  Mach  Numbers  ,  AIAA  J., 
31(1993),  pp.  49-56. 


18 


log(res) 


2.00  n 

j 

d 


0.00 


Hgure  1:  Residual  for  the  moineutum  equation.  Flow  about  a  NACA0012, 
Moo  =  0.03, a  =  4®  with  a  98*22  C  grid.  Graph  (1)  is  the  preconditioned 
solution  and  (2)  is  without  preconditioning. 
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figure  2;  Same  as  figure  I  but  with  Mc^  -  0.80,  a  1-25 
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